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Abstract 

Simulation-based verification algorithms can provide formal safety guarantees for nonlinear 
and hybrid systems. The previous algorithms rely on user provided model annotations called dis¬ 
crepancy function, which are crucial for computing reachtubes from simulations. In this paper, 
we eliminate this requirement by presenting an algorithm for computing piece-wise exponen¬ 
tial discrepancy functions. The algorithm relies on computing local convergence or divergence 
rates of trajectories along a simulation using a coarse over-approximation of the reach set and 
bounding the maximal eigenvalue of the Jacobian over this over-approximation. The resulting 
discrepancy function preserves the soundness and the relative completeness of the verification 
algorithm. We also provide a coordinate transformation method to improve the local estimates 
for the convergence or divergence rates in practical examples. We extend the method to get 
the input-to-state discrepancy of nonlinear dynamical systems which can be used for composi¬ 
tional analysis. Our experiments show that the approach is effective in terms of running time 
for several benchmark problems, scales reasonably to larger dimensional systems, and compares 
favorably with respect to available tools for nonlinear models. 


1 Introduction 

Verifying and falsifying nonlinear, switched, and hybrid system models using numerical simulations 
have been studied in detail [10, 17, 4, 14, 9]. The bounded time safety verification problem for a 
given model is parameterized by a time bound, a set of initial states, and a set of unsafe states and 
it requires one to decide if there exists a behavior of the model that reaches any unsafe set from 
any initial state. The simulation-based procedure for this problem first generates a set of numerical 
approximations of the behaviors from a finite sampling of the initial states. Next, by bloating these 
simulations by an appropriately large factor it computes an over-approximation of the reachable 
states from the initial set. If this over-approximation proves safety or produces a counter-example, 
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1054247 and NSF CSR 1016791) and the Air Force Office of Scientific Research (AFOSR YIP FA9550-12-1-0336). 
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then the algorithm decides, otherwise, it draws more samples of initial states and repeats the earlier 
steps to compute more precise over-approximation. With post-processing of the reachtube over¬ 
approximations this basic procedure can be utilized to verify termporal precedence [11] and richer 
classes of properties [7]. 

In order to make this type of procedure sound, the bloating factor should be chosen to be large. 
Specifically, it should be large enough to make each bloated simulation an over approximation 
of the reachable states of the system not only from the sampled initial state, but also from a 
large enough neighborhood of that state so that the union of these neighborhoods cover the entire 
set of initial states. On the other hand, to make the procedure complete, or at least relatively 
complete modulo the precision of the machine, it should be possible to make the error due to 
bloating arbitrarily small for any point in time. These two opposing requirements are captured in 
the definition of a discrepancy function of [10]: For an n-dimensional dynamical system, it is any 
function fj : x M>o —M>o, such that (a) it gives an upper-bound on the distance between any 

two trajectories and of the system —(,{x',t)\ < j3{x,x',t), and (b) it vanishes 

as x approaches x'. Simply using the Lipschitz constant of the dynamic function gives one such 
bound, but it grows exponentially with time even for some incrementally stable models [2]. 

In [10], it is observed that the notion of a contraction metric [19] gives a much tighter bound 
and it provided heuristics for finding them for some classes of polynomial systems. Sensitivity 
analysis approach gives strict error bounds for linear systems [9], but for nonlinear models the 
bounds are less clear. We present a more detailed overview of related work in Section 6. This paper 
fills this gap by providing a subroutine that computes a local version of the discrepancy function 
which turns out to be adequate and effective for sound and relatively complete simulation-based 
verification. This subroutine, ComputeLDF, itself uses a Lipschitz constant and the Jacobian of the 
dynamic function (the right hand side of the differential equation) and simulations of the system. 
The Lispchitz constant is used to construct a coarse, one-step over-approximation of the reach set 
of the system along a simulation. Then it computes an upper bound on the maximum eigenvalue 
of the symmetric part of the Jacobian over this over approximation, using a theorem from matrix 
perturbation theory. This gives an exponential bound on the distance between two trajectories, 
but roughly, the exponent is the best it can be as it is close to the maximum eigenvalue of the 
linear approximation of the system in the neighborhood. 

We propose two practical extensions of this approach. First, we show that a linear coordinate 
transformation can bring about exponential improvements in the estimated distance. Secondly, 
we propose a technique for computing input-to-state discrepancy functions for analyzing composed 
systems and systems with bounded nondeterministic inputs. We report the results from a number 
of experiments performed with a prototype implementation of this approach applied to safety 
verification. 

2 Background 

2.1 Notations 

For a vector x E M”, ||x|| is the /^-norm of x and Xi denotes its component. For J > 0, 
Bs{x) = {x' E M"" I |]x' — x|| < J}. For a set 5 C M"-, Bs{S) = \Jx^sB 5 {x). Let S®Bi{Q) represents 
the Minkowski sum of S and B^{t)). Therefore, S 0 1?5(0) = Bs{S). For sets 81,82 C M"’, 
hull{ 8 i, 82 ) is their convex hull. The diameter of a compact set 8 is dia{ 8 ) = sup^^.^ X 2 £S ll^i “® 2 ||- 
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A continuous function / : MT' —>■ M is smooth if all its higher derivatives and partial derivatives 
exist and are also continuous. It has a Lipschitz constant L > 0 if for every x,x' G M”, \\f{x) — 
f{x')\\ < L\\x — x'||. A function / : M>o —5- M>o is a class K, function if it is continuous, strictly 
increasing, and /(O) = 0. 

We denote the transpose of a matrix A by The conjugated transpose of A is the matrix A^ 
obtained by replacing each entry in with its complex conjugate. 

Given a differentiable vector-valued function / : x —>• W^, the Jacobian Jf of / is the 

matrix-valued function of all the first-order partial derivatives of /. Let fi,i = ^ ■ ■ ■ n : M>o 

be the scalar components of /. The Jacobian of / is: {Jf{x))ij = The symmetric part of the 

Jacobian of f matrix is defined as ^(J/(x) -|- Jf'^{x)). 

For an n X n matrix A, ||A|| represents the Z^-norm of A: ||A|| = \/ Xmax{A^A). If Vx G 
x^Ax < 0, then we say A is negative-semidefinite, and write A ^ 0. We write A < B if A — B < 0. 

2.2 Safety Verification Problem 

Consider an n-dimensional autonomous dynamical system: 

x = fix), (1) 

where / : M”’ ^ M"' is a Lipschitz continuous function. A solution or a trajectory of the system is a 
function ^ : M"’ x M>o —>■ M"' such that for any initial point xq G M” and at any time t > 0, ^(xo,t) 
satisfies the differential equation (1). 

The bounded-time safety verification problem is parameterized by: (a) an n-dimensional dynam¬ 
ical system, that is, the function / defining the right hand side of its differential equation, (b) a 
compact set 0 C M” of initial states, (c) an open set U C M** of unsafe states, and (d) a time 
bound T > 0. A state x in M"' is reachable from 0 within a time interval \ti,t‘2\ if there exists an 
initial state xq G 0 and a time t G [ti,t 2 ] such that x = ^{xo,t)- The set of all reachable states 
in the interval [ti,t 2 ] is denoted by Reach(0, [fi,f 2 ]). If = 0 then we write Reach(t 2 ) when set 
0 is clear from the context. Given a bounded-time safety verihcation problem, we would like to 
design algorithms for deciding if any reachable state is safe, that is, if Reach (T) n U = 0. If there 
exists some e > 0 such that JIe(Reach(T)) nU = 0, we say the system is robustly safe. A sequence 
of papers [10, 11, 9] presented algorithms for solving this problem for a broad class of nonlinear 
dynamical, switched, and hybrid systems. In the remainder of this section, we present an overview 
of this approach. (Figure 1). 

2.3 Simulations, Reachtubes and Annotations 

The algorithm uses simulation oracles that give sampled numerical simulations of the system from 
individual initial states. 

Definition 2.1. A {xo,t,€,T)- simulation of the system described in Equation (1) is a sequence of 
time-stamped sets (iZoAo); (RiAi) • • • > {RnAn) satisfying: 

(1) Each Ri is a compact set in with dia{Ri) < e. 

(2) The last time tn = T and for each i, 0 < fj — ti-i < r, where the parameter r is called the 
sampling period. 
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(3) For each ti, the trajectory from xq at U is in Ri, i.e., G Ri, and for any t G 

the solntion ^{xo,t) G hull{Ri-i, Ri). 

Simnlation engines generate a sequence of states and error bounds using numerical integration. 
Libraries like CAPD [5] and VNODE-LP [20] compute such simulations for a wide range of nonlinear 
dynamical system models and the Ri’s are represented by some data structure like hyperrectangles. 

Closely related to simulations are reachtubes. For a set of states D C M"-, a {D,T,T)-reachtube 
of (1) is a sequence of time-stamped sets {Rq, 0), {Ri,ti) ..., (i?„, satisfying: 

(1) Each Ri C M” is a compact set of states. 

(2) The last time tn = T and for each i, 0 <ti — U-i < r. 

(3) Eor any xq G D, and any time t G the solution ^{xo,t) G Ri- 

A reachtube is analogous to a simulation from a set of states, but they are much harder to compute. 
In fact, an algorithm for computing exact reachtubes readily solves the safety verification problem. 

The algorithms in [10, 16] require the user to decorate the model with annotations called 
discrepancy functions for computing reachtubes. 

Definition 2.2. A continuous function (3 : x M” x M>o —M>o is a discrepancy function of the 

system in Equation (1) if 

(1) for any pair of states x, x' G M”, and any time t > 0, 

||^(x,t) - ^(x',t)ll < P{x,x',t),and (2) 

(2) for any t, as x ^ x', /3(.,., f) ^ 0, 

If the function fi meets the two conditions for any pair of states x,x' in a compact set K then 
it is called a K-local discrepancy function. 

The annotation (3 gives an upper bound on the distance between two neighboring trajectories 
as a function of their initial states and time. Unlike incremental stability conditions [2], the second 
condition on (3 does not require the trajectories to converge as time goes to infinity, but only as the 
initial states converge. Obviously, if the function / has a Lipschitz constant L, then I3{x,x',t) = 
jjx — meets the above criteria. In [10, 16] other heuristics have been proposed for finding 

discrepancy functions. As will be clear from the following discussion, the quality of the discrepancy 
function strongly influences the performance of the simulation-based verification algorithm. [10, 
16, 17] need user provided discrepancy and simulation engines to give verification of bounded time 
safety and temporal precedence properties. In this paper, we will present approaches for computing 
local discrepancy functions that unburdens the user from finding these annotations. 

2.4 Verification Algorithm 

The simulation-based verification algorithm is shown in Eigure 1. It takes as input some finite 
description of the parameters of a safety verification problem, namely, the function /, the initial set 
0, the unsafe set U, and the time bound T. It has two main data stuctures: The first, C returned 
by function Partition, is a collection of triples {0,6, e) such that the union of all the h-balls around 
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1: Input:©, U, T 

2: (5 •(— dia{Q)]e ■(— eo;C •(— 0,7^ •(— 0; //eo is a small constant 
3: C •(— {Partition{Q,6),6,€) 

4: while C / 0 do 
5: for {9, 6,e) ^ C do 

6: "0 •<— Simulate{9, r, e, T) 

7: /3 •(— ComputeLDF{'ijj,Jf,Lf,5,e) 

8: D ^ V'®/3 

9: if Z7 n U = 0 then 

10: C^C\{(0,(5,e)};7^^7^UZ) 

11: else if 3k, Rk C U then 

12: return (UNSAFE, 7^) 

13: else 

14: C^C\{(0,5,e)} 

15: C C U Partition{Q n Bg{6), (^,..., ^), |) 

16: end if 

17: end for 

18: end while 

19: return (SAFE, 7^) 


Figure 1: Verification Algorithm 


the 0’s completely cover the initial set 0. The second data structure TZ incrementally gets the 
bounded-time reachtube from 0. 

Initially, C has a singleton cover {6o,6o,eo) such that (5o = dia{Q), 0 C Bs(,{6o), and cq is a 
small constant for simulation precision. 

In the while-loop, this verification algorithm iteratively refines the cover of 0 and for each 
(0, 6, e) in C, computes over-approximations of the reachtube from Bs{9). The higher-level structure 
of the algorithm is familiar: if the reachtube from B^{6) proves to be safe, i.e., disjoint from U, 
then the corresponding triple is removed from C (Line 10). If part of the reachtube from Bs{6) 
overlaps with U, then the system is declared to be unsafe (Line 12). ©therwise, a finer cover of 
Bs{9) is created, and the corresponding triples with finer parameters are added to C. 

Here we discuss the reachtubes computed from discrepancy and simulations. For each (0, 6, e) 
in C, a (0, r, e, r)-simulation ijj, which is a sequence of {(77*, fj)}, is generated. Note that -ip contains 
the trajectory from 0, ^(0,7), t G [0,T]. Then we bloat each 77* by some factor (Line 8) such 
that the resulting sequence contains the reachtube from Bs{6). It is shown that this bloated 
simulation is guaranteed to be an over-approximation of Reach (77^(0), T) and the union of these 
bloated simulations is an over-approximation of Reach(0,r). Therefore, the algorithm is sound. 
Furthermore, the second property of (5 ensures that the reach set over-approximations become 
tighter and tighter as we make 5 smaller and smaller. Finally it will return “SAFE” for robustly 
safe reachtubes or find a counter example and return “UNSAFE”. For user defined discrepancy 
function, the factor is obtained by maximizing /3(0,0,7) over 0 G 77^(0) and t G [7j_i,7j]. 

Indeed this is the approach taken in the algorithm presented in [10]. In this paper, we will ana¬ 
lyze in detail the ComputeLDF subroutine which computes a local version of discrepancy function 
automatically. 
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The following results from [10] state two key properties of the algorithm. Although in [10] /3 
was defined globally, it is easy to check that the local version still satisfies them. 

Theorem 2.3. The Algorithm in Fig.l is sound, that is, if it returns “SAFE” then the system is 
safe; when it returns “UNSAFE” there exists at least one execution from 0 that is unsafe. The 
Algorithm is relatively complete, that is, if the system is robustly safe, the algorithm will terminate 
and return “SAFE”. If any executions from 0 is unsafe, it will terminate and return “UNSAFE”. 

3 Local discrepancy function 

In this section, we present the analysis of ComputeLDE algorithm. This algorithm computes 
a special type of local discrepancy in terms of time-varying exponential functions that bound 
from above the distance between two trajectories starting from a compact neighborhood. Roughly 
speaking, it computes the rate of trajectory convergence or divergence for an interval of time 
instances. 

Definition 3.1. Consider a compact set C C M"' and a sequence of time points 0 = to < < 

t 2 < ... < tk = T. For Vxi,X 2 G C,Vt G [0,^], a piece-wise exponential discrepancy function 
/3 : C X C X [0,T] ^ R>o is defined as 
(3 {xi,X2,t) = 

{ \\xi-X 2 \\, At = to, 

if t G iti-i,ti], 

where 5[1],..., b[k] are real constants. 

From the definition, we can immediately get that / 3 {xi,X 2 ,t) = \ |xi—X 2 I 
i = 1,..., k, where fj-i is the largest time point in the sequence before t. 

3.1 ComputeLDF Algorithm 

Figure 2 shows the pseudocode for ComputeLDF used in Line 7 of the verification algorithm. 
ComputeLDF takes as input a parameter 5, an error bound for simulation e, the Lipschitz constant 
Lf, the Jacobian matrix Jf of function /, and a (0, r, e, T)-simulation = {(R*, tj)}, f = 0,1,... , A:. 
It computes a piece-wise exponential local discrepancy function (LDF) for the compact set Bs{Ro) 
and for the time points tQ,... ,tk. and returns it as an array of exponential coefficients b. 

The algorithm starts with the initial set Bs{Ro) and with A = 5. In each iteration of the for- 
loop it computes exponent b[i] corresponding to the time interval [ti-i,ti]. In the iteration, A 
is updated so that RaC-Rj-i) is an over-approximation of the reachable states from Bs{Ro) at 
(Lemma 3.8). In Lines 8 and 9, a set S is computed by bloating the convex hull hull{Ri-i, Ri) by a 
factor of d = The set S will later be proved to be a (coarse) over-approximation 

of the reachtube from B^{Ri-i) over the time interval [L_i,ti] (Lemma 3.2). In Lines 10-13 an 
upper bound on the maximum eigenvalue of the symmetric part of the Jacobian over the set S, is 
computed as b[i] (Lemma 3.6). Then A is updated as (A -|- e)e^W(L-L-i) next iteration. 


6 


1: Input: 

2: A ■(— 6,b •<— zeros(k) 

3: for i = l:k do 
4: T ^ti- tj_i 

5: d ^ (A + e)e^f^ 

6: S •(— hull{Ri-i, Ri) © i?d(0) 

7: J ■(— Jf {center{S)) 

8: A ma.x{eig{J + J^)/2) 

9: errors xes\\{Jf{x) + Jj {x)) - {J + J^)\\ 

10: b[i] A + error/2 

11: A ^ (A + e)e'’W^ 

12: end for 
13: return b 


Figure 2: Algorithm ComputeLDF. 

3.2 Analysis of ComputeLDF 

In this section, we will prove that ComputeLDF{'ilj,Jf,Lj-, 6,e) returns a piece-wise exponential 
LDF of the system in Equation (1), for the compact neighborhood Bs{Ro), and the sequence of the 
time points in the simulation ■0. We establish some lemmas to prove the main theorem. First, we 
show in Lemma 3.2 that in the iteration of the loop, the computed S is an over-approximation of 
the set of states that can be reached by the system from Bj\{Ri-i) over the time interval 

Lemma 3.2. In iteration of the loop of ComputeLDF, Reach {B^{Ri-i), [ti-i,ti\) C S. 

Proof. Let f{9,t) denote the actual trajectory from 9, where 9 is the initial state of By Defini¬ 
tion 2.1 for ijj, it is known that 9 £ Rq and Vi = 1, ... ,k,f,{9,ti) E i?j. 

For a fixed iteration number i, consider state x = f^{9,ti-i) E Ri-i from Definition 2.1. We 
know that for any t E E hull{Ri-i, Ri). Now consider another state x' E 

Since Lj is the Lipschitz constant of /, using Gronwall’s inequality we have that ||^(x, t)—f,{x',t)\\ < 
||x—x'||e'^7V-L-i)^ Since ||x—x'|| < A+e, ||^(x,t)—^(x',t)|| < (A+e)e-^7V-L-i)_ Therefore, ^(x', t) E 
hull{Ri-i, Ri) © (0) = S. Because x' is arbitrarily selected from the 

lemma is proved. ■ 

Next, using the generalized mean value theorem (Lemma 3.3), we get that in the iteration, 
the computed b[i] in Line 13 is the exponential divergence (if positive) or convergence (negative) 
rate of the distance between any two trajectories starting from B^{Ri-i) over time [ti-i,ti]. 

Lemma 3.3. For any continuously differentiable vector-valued function / : ^ and x,r E 

R", 

f{x + r)-f{x)= Jf{x + sr)ds^ -r, (3) 

where the integral is component-wise. 

Next, we will use a well-known theorem that gives bounds on eigenvalues of perturbed symmetric 
matrices, the proof of which uses the Courant-Fischer minimax theorem. The complete proofs of 
Lemma 3.3 and Theorem 3.4 can be found in the appendix. 


7 



Theorem 3.4. If A and E are n x n symmetric matrices, then 


K{E) < \k{A + E)- Xk{A) < Ai(E), 

where Aj(-) is the largest eigenvalue of a matrix. 

Corollary 3.5. If A and E are n x n symmetric matrices, then 

\Xk{A + E)-Xk{A)\<\\E\\. (4) 

Since A is symmetric, ||j 4|| = Ajnax(^^^) = max(|A(T)|). From Theorem 3.4, we have 
|Afc(^ + £') — Xk{A)\ < max{ I An (F^) I, |Ai(£^)|} = ||F^||. If E{x) is a matrix-valued function: R” —)■ 
]^nxn j^^g^pg g stutc X G R”’ to a matrix E{x), and every component of E(x),eij{x) : R”’ ^ R is 
continuous over some compact closed set S, then we can ||Fl(x)|| over S by each term eij(x), 
\eij{x)\ over S. Let x&s{\eij{x)\) be denoted by Cij, then we know , ||Fl(x)|| < 

\JY11=i Using Corollary 3.5, we next show in Lemma 3.6 that h[i\ calculated in Line 13 

bounds the eigenvalues of symmetric part of Jacobian matrix over S. 

Lemma 3.6. In the i^^ iteration, for Vx G 5* : Jj{x) + Jf{x) < 2h[i]I. 

Proof. Let S be the set computed in Line 9 and J be the Jacobian evaluated at the center sq of S. 
Consider any point x G 5. We define the perturbation matrix E{x) = Jj{x) + Jf{x) — (J^ J- J). 
Since Jj (x) -I- J/(x) and + J are symmetric matrices. Corollary 3.5 implies that Xmax{Jj(x) + 
Jf(x)) — Xmax{J'^ + F) < ||Fl(x)||. The error term computed in Line 12 is the upperbound on 
||Fl(x)||. Therefore, Xmax{Jj(x) + Jf{x)) < Xmax{J"^ + J) + error. In Line 13 set b[i] equals to 
^max{{J'^ + F)/2) J- error/2. Thus, XmaxiJj (x) J- Jf{x)) < 26[i], which immediately indicates that 
'Ix £ S : jJ (x) J- J/(x) R 2b[i]I. ■ 

By Lemma 3.3 and Lemma 3.6, we can prove as in Lemma 3.7 that b[i] calculated in Line 13 is 
the exponential rate of divergence or convergence of two trajectories starting from B,\{Ri-i) over 
the interval [L_i,fi]. 

Lemma 3.7. In the i^^ iteration, for any two states xi,X 2 G at time tj-i, and any time 

t£ [ti-i,ti], ||^(xi,t) - ^(X 2 ,t)|| < ||xi - X 2 ||e'’WF-L-i)_ 

Proof. Let us fix the iteration i and two states xi,X 2 G Bj\{Ri-i). From Lemma 3.2 it’s can be 
seen that for any t G L], ^(xi,t) G S,(^{x 2 ,t) G S. Define y{t) = ^(x 2 , f) — ^(xi, f). For a fixed 
time t, from Lemma 3.3 we have 

y{t) = iix 2 ,t) - i{xi,t) = f{f{x 2 ,t)) - f{f{xi,t)) 

= Jf{i{xut) +sy{t))d^y{f). (5) 

Since S is the Minkowski sum of two convex sets hull{Ri-i, Ri) and (0), it is also 

convex. Recall that ^(xi,f), f{x 2 ,t) G S, and for any s G [0,1], f{xi,t) + sy{t) C S. 

Differentiating ||y(f)|p, we have 




dimr 

dt 


= iF {t)y{t) + y'^ {t)y{t) 

= y^{t) Jf + sy{t))ds^ y{t) 

+ y^it)( [ Jf{^ixi,t) + sy{t))ds]y{t). 


Using Lemma 3.6, we know 

Vx E S, Jj{x) + J/(x) :< 2b[i]I. 

Thus, we can bound (6) 

= ‘^b[i]y'^{t)y{t) 

= ni\\\yit)f- 

Integrating both sides over U-i to any t E we have 

ln(|| 2 /(t)f)-ln(|| 2 /(U_i)f)< 26 [i](t-U_i) 

- ^(a^2,i)|| < ||xi - 


(6) 


(7) 


Up to this point all the lemmas were statements about a single iteration of the for -loop, next 
we show that in iteration of the loop, used in Lemma 3.2 and 3.7 is the reach set from 

Bs{Ro) at time U. 

Lemma 3.8. For Vz = Reach [U,U]) U BAi{Ri), and Reach U 

hull{Ri-i, Ri) 0 i?A'(0)) where Aj is A after Line 14 is executed in the z*^ iteration, and A'- = 
max{Ai, Ai_i 0 e} 

Proof. In this proof, let denote the trajectory from 6 . From the Definition 2.1 for ip, we 

know that 6 £ Rq and Vz = 1,..., A:, ^(0, ti) E Ri- Let Si denote S after Line 9 is executed in the z*^ 
iteration. The lemma is proved by induction on z. Note that the initial set is Bs{Rq), and before 
the for-loop, Aq is set as <5. 

When z = 1, we already have Reach (i? 5 (i?o)) [Aoj^o]) = Bs{Ro) = 7?Ao(7?o)- 

Lemma 3.2 indicates that Vt E [to,U], Reach (i?Ao(7?o)) [Ao,Ai]) U S. And consider state x = 

9 E Rq, we also know f,{x,t) E hull{Ro, Ri) and ^(x,ti) E Ri. From Lemma 3.7, it follows that for 
Vx' E BAo{Ro),'^t E [to,U], 

||^(x,t) — ^(x',t)|| < ||x — x'||e^Wd-*o), 
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And at Line 14, Ai •(— (Ao+e)e^I^l(*i Since 6[1] could be positive or negative, ||x — 

max{Ai, Aq + e}. Therefore, 

Reach (^5 (i^o), [ioTl]) ^ ^l) ® 5max{Ai,Ao+e}(0)> 

and at time ti, is at most Ai distance to ^(x,ti) G Ri, so Reach (i?5(i2o)) [^iTi]) = 

Reach(5 ao(-R o), [iiTi]) ® Bai{Ri). 

Assuming that the lemma holds for i = m—1, we have Reach ® -Ra^_i (-Rm-i)- 

Next we prove the lemma holds for z = m as well. Consider state x = G Rm-i, 

Vt G by definition it follows that G hull{Rm-i, Rm) and ^{x,tm) G Rm- Vx' G 

^A^_i(i?m-i),Vt G from Lemma 3.7 

||^(x,t)-e(x',t)|| < ||x-x'||e^Hh-t—i). 

Note at Line 14, A^ <- (A^-i + Therefore, Reach(S a„_i), [im-iTm]) ® 

hull{Rm-i, Rm) ® -Bmax{A^,A^_i+e}(0). And at time tm, i{x',tm) is at most Am distance to 
^{Xitm) S Rm- Hence, Reach(HA,^_j(7?m—i)) [^mTm]) ® RAmi^m)- Recall that Reach(^^(R q), [tm—iTm—i]) 
® RAm—i^Bm—l)i thus Reach (Rq) , [tmTm]) ® HA^(Rm)' ® 

U^^^{/ittR(Rj_i, Rj) © Ra'(O)} contains the (R 5 (Ro), r, r)-reachtube of the system. Line 8 of 
the algorithm in Figure 1 is computed in this way. Now we are ready to prove the main theorem. 

Theorem 3.9. The items in array b computed by 

ComputeLDF are the coefficients of a Bs{Rq) -local piece-wise exponential discrepancy function 
(Definition 3.1). 

Proof. First of all consider any time t G [toTi] and any two states: xi, X 2 G Bs{Rq). By Lemma 3.7, 

||^(xi, t) —^(x 2 , t)|| < jjxi — X 2 ||e^[^]*'^“^°). Then consider t G [ti,t 2 ]- By Lemma 3.8 we know at time 
ti, f,{xi,ti) and (,{x 2 ,ti) are all contained in Rai(Ri)) so we can use Lemma 3.7 such that for any 
timet G [ti,t 2 ], ||C(a^iT)-?(a;2T)ll < M{xi,ti)-f,{x 2 ,ti)\\e^^‘^^'^^~^^'> < jjxi 

The procedure above can be performed iteratively as follows. For any time t G [ti_i,ti], by 
lemma 3.8 we know at time tj_i, ^(xi,tj_i) and ^(x 2 ,tj_i) are all contained in RAi_iBy 
Lemma 3.7 it follows that 

||^(xi,t)-e(x2,t)|| < ||^(xi,L_i)-C(x2,t,_i)||e'W(‘-‘»-i) 

Next we will prove that 

/3(xi,X2,t) = jjxi — is a valid LDF. 

In Lines 10-13, because J is a real matrix, the maximum eigenvalue A of (J^ +J)/2 is bounded. 

Assume that each component of E{x) = Jj (x) + J/(x) — — J is continuous over the closed set S, 

, so the “error” term is also bounded. Therefore, each 6[z] is bounded. So Vt G [ti_i, tj], t = 1,..., A;, 

3N < oo, such that is bounded by N from the above. 

As xi —>■ X 2 , obviously, 

\\xi - ^ 0. 
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And for any e > 0, 35 = e/N > 0, such that Vxi,X 2 G Bs{Ro) and ||xi — X 2 II < 5, it follows 

||xi - <e/N ■ N = e. 

So l3{xi,X2,t) = ||xi — is a i? 5 (i?o)-local piece-wise discrepancy func¬ 
tion and the array b contains the corresponding coefficients. ■ 


3.3 Coordinate transformation 

In this section, we will discuss the issue that the upper bound of the symmetric part of the Jacobian 
computed in Lines 10-13 may introduce loss in precision. We propose a a strategy to reduce this 
loss by first performing a coordinate transformation. Consider a simple linear system: 


X = 


0 3 

-1 0 


X, 


( 8 ) 


which has eigenvalues zb\/3i and thus its trajectories oscillate. The symmetric part the of the 
0 1 


Jacobian is 


1 0 


with eigenvalues ±1, which gives the exponentially growing discrepancy 


with 6 = 1. In what follows, we will see that a tighter bound can be obtained by first taking 
linear transformation of x. The following is a coordinate transformed version of Lemma 3.7. The 
coordinate transformation matrix P can be any n x n real invertible matrix, and the condition 
number of P is ||P|| ||P“^ ||. 

Lemma 3.10. In iteration of the loop, for any xi,X2 G PA(Pi-i), and any t G 

||^(xi,t) - I(x 2 ,t)|| < K\\xi - 

where Xmax{S) is the upper bound of ^(J/ (x) -|- Jf{x)) over the set S, Jf{x) = PJf{x)P~^, and 
K is the condition number of P. 

Proof. Let z{t) = f^{x2,t) — ^(xi,t) and y{t) = Pz{t) From (5) get: 

y{t) = Pz{t) 

= P 

= P 


Jf{f,{xi,t) -I- sz{t))ds z{t) 


Jf{f{xi,t) + sz{t))ds \ P ^y{t) 


+ sz{t))dsj y{t). 

Since for all x G S,Jf (x) + Jf{x) P Amax(*S')/ and Vs G [0, l],^(xi,f) -|- sz{t) C S, we have 

d\\y{t)r 


dt 


<2XmUS)\\y{tW, 
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which leads to: Vt G 


||y(t)|| < ||y(t,_i)||e^—(9) 
Substituting (9) in z{t) = P~^y(t): 

||^(t)|| < ||P-1l|yWII 

= cond(P)||2;(ti_i)||e^’"“''^'^^*'*“*‘“^^. (10) 


This shows that the distance can be bounded in the same way for the transformed system with 
a (possibly much smaller ) Xmax{S) but with an additional multiplicative cost of cond{P). 

Let J = PJP~^ the real Jordan form which looks like: 


Ai 

e 

0 

0 

0 

0 

Ai 

0 

0 

0 

0 

0 

A 2 

0 

0 

0 

0 

0 

A 3 

c 

0 

0 

0 

—c 

A 3 


where 2e < Ai and Ai, A 2 , A 3 ± ci are the eigenvalues of J. There could be several more blocks like 


Ai e 
0 Ai 


,A 2 and 


A 3 c 
-c A 3 


in general. We use the matrix P as the coordinate transformation 


matrix for Jf{x). In this approach the eigenvalues of ^( J + J^) are Ai + A 2 , A 3 , which preserve 

the original eigenvalues to some extent. 

In the previous example ( 8 ), the Jacobian matrix is constant, and the discrepancy function 
without coordinate transformation is: 


||^(xi,t) -^(X2,t)|| < ||xi -X2||e* 


1 

3 


0 

v^’ 

_-V3 

V3_ 

as the coordinate transformation matrix, J = PJP ^ = 

. 

0 


and the discrepancy function with coordinate transformation is 

||^(xi,t) -^(X2,t)|| < V3||xi -X2 


In practice, the coordinate transformation can be made for longer time interval where 

A: > 2, to reduce the multiplicative error term ])([ cond{P[i]). 
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1 : Input: ip,Jf ,Lf ,6, €,step 

2: A ■(— 6 ,b •(— zeros(k), K •(— zeros(ceil(k/step)) 

3: for j = l:step:k do 

4 ; [V, D] = JordanDecom(average(i?j,...,i?j+step-i)) 

5: K{ceil{k/step))) ^ cond{V) 

6: for i = j:j+step-l do 

7: T ^ti- ti-i 

8: d ^ (A + e)e^7'^ 

9 : S hull{Ri-i, Ri) (B B(i{0) 

10 : J ^ VJf {center{S))V-^ 

11: A •(—max(ei(7(J + j'^)/2) 

12: error upper xesW^((.Jf i^) + i^)) “ 

13: b[i] ^ A + error /2 

14: A ^ (A + 

15: end for 

16: A K{ceil{k / step)))A 

17: end for 
18: return b, K 

Figure 3: Algorithm ComputeLDF to coordinate transformation. 

4 Local Input-State Discrepancy 

Large and complex models of dynamical system are created by composing smaller modules or sub¬ 
systems. Consider a dynamical system A consisting of several interacting subsystems Ai,..., A^r, 
that is, the input signals of a subsystem A* are driven by the outputs (or states) of some another 
component Aj. Let’s say that each Aj is n-dimensional which makes A nA-dimensional. One 
way of achieving scalable verification of A is to exploit this compositional structure and somehow 
analyze the component A/s to infer properties of A. 

In [17], the notion of input-to-state (IS) discrepancy was introduced to address the problem of 
finding annotations for large models. It is shown that if we can find input-to-state (IS) discrepancy 
functions for the individual component Aj, then we can construct a reduced A-dimensional model 
M such that the executions of M serve as the discrepancy of the overall system. Thus, from IS- 
discrepancy for the smaller Aj models and simulations of the A-dimensional system M, we are 
able to verify A. This has the beneficial side-effect that if the Aj’s are rewired in a new topology, 
then only the reduced model changes [16]. However, [17] still assumes that the user provides the 
IS-discrepancy for the smaller modules. In this section, we will show the approach used in previous 
section can be used to get IS discrepancy function for Lipschitz continuous nonlinear subsystems 
Aj. Furthermore, it gives an over-approximation of the reachsets with nondeterministic bounded 
inputs. 

4.1 Defining Local IS Discrepancy 

Consider a dynamical system with inputs: 

x = f{x,u) ( 11 ) 
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where / : M"" x —)• M"" is Lipschitz continuous. For a given input signal which is a integrable 
function v : [ 0 ,oo) ^ and an initial state xq € M”, a solution (or trajectory) of the system is a 
function ^ x M>o ^ M”' such that ^(xq, 0 ) = xq and for any time t > 0 , i{x, t) = /(^(x, t),v{t)). 

First, we give the original definition of IS discrepancy function for the system in (11). Here U 
is the set {u\u : [ 0 ,oo) —>■ of all input signals. 

Definition 4.1. A pair of uniformly continuous functions j3 : M>o xR>o —)■ M>o and 7 : M>o —>• R>o 
is called C-local input-to-state discrepancy if 

(1) (3 is of class K, with respect to its first argument and 7 is also of class /C, 

(2) for any pair of initial states x, x' G C, any pair of input signals u, u' G U, and t G M>o: 

||C(x,t) - C(x',t)|| </3(||x - x'||,t) + / 7 (||u(s) - u'(s)||)ds. ( 12 ) 

Jo 

For a bounded, compact set I C R^. A family of bounded time input signals over I is the 
set U{T) = {u\u : [ 0 ,r) —>■ X} of integrable functions. We denote Reach(AT,W(X), [ti,t 2 ]) as the 
reachable states of the system from compact set K with input set U{X) over [ti,t 2 ]- Next, we 
introduce an inductive definition of IS discrepancy for inputs over compact neighborhoods. 

Definition 4.2. Consider compact sets K G R”’,X G R^ and a sequence of time points 0 = to < 
ti < t 2 < ... < tk = T. For any pair of initial states xi,X 2 G K, any pair of input signals 
ui,U 2 G U{X), the {K,U{X))-local IS discrepancy function a : x UiX)"^ x R>o —>■ R>o is defined 

as: 

a(xi,X 2 ,ttl,tt 2 ,t) = 

||xi - X 2 II, if t = toj 
< q;(xi, X 2 , ui, n 2 , 

^ ll^ii(T) - U2{T)\\dT,\i t G (ti-i,ti], 

where a[l],... , a[k\, Af[l],... , M\k] are real constants. 

4.2 Algorithm for Local IS Discrepancy 

The approach to find [K^IA{X))-\oca\. IS discrepancy function is similar to ComputeLDF algorithm, 
which also uses a for -loop to compute the coefficients a[i\ and M\i]. The only changes are 1) in 
Line 9 S should be computed as in Lemma 4.3, 2) in Line 14 A should be updated as in Lemma 
4.5. Next we illustrate this process in more detail. First, we use Lipschitz constant to get a coarse 
over-approximation of Reach(iL, Z//(X), [tj_i,L]) parallel to Lemma 3.2. Let I = dia{X). 

Lemma 4.3. In iteration of the for -loop, Reach{BA{Ri-i), U{X),[ti-i,ti]) C S, where S = 
hull{Ri-i, Ri) 0 BA'iRi) and A' = (A 0 0 ILfC^f^^n, n = U — A_i. 

Two trajectories starting from xi,X 2 G R”' at A_i, with ui,U 2 G U{X) as inputs respectively, 
their distance at time t, ||^(xi, t)—^(x 2 , t)|| < ||xi—X 2 ||(e'^J'(*“**-i ))0 ^ ||ui(r) — U 2 (T)||(ir. 

The lemma directly follows this inequality. 

Next we give a one step IS discrepancy function in Lemma 4.5. Before proving it, we need 
another generalized form of mean value theorem: 
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Lemma 4.4. For any continuous and differentiable function / ; M"' x ^ M”, f(x + r,u + w) — 
f{x,u) = 

(fo + sr,u + w)ds^ r+ ^ Ju(x, u + Tw)dT^ w, where Jx = and Ju = are the 

Jacobian matrices of / with respect to x and u. 

Proof. The lemma follows by writing f{x + r,u + w) — f{x, u) = f{x + r,u + w) — f{x, u + w) + 
f{x, u + w) — f{x, u) and then invoking Lemma 3.3. ■ 


Lemma 4.5. Consider the iteration of the loop for a dynamic system (11). Let x, x' G 
and (,{x,t), be the trajectories starting from x and x' with input ui{t),U2{t) G respec¬ 
tively, where t G Then, 

||C(x,t)-e(x',t)|| < ||x - 

-|- f \\ui{t) — U2 {T)\\dT, (13) 

J ti — \ 


where a = Xmax{S) + Xmax{S) is the upperbound of the eigenvalues of the symmetric part of Jx 
over S', and M = max (|| Ju(^(x,t),u)||). 

u<3J{X) 

Proof, let y{t) = f^{x',t) — f,{x,t) and v{t) = U2{t) — uift). For a fixed time t, using Lemma 4.4 


y{t) = C{x',t) - f{x,t) 

= f{Cix',t),U2{t)) - f{f,{x,t),Ul{t)) 


Jx{^ix,t) + sy{t),U2{t))ds^ y{t) 


+ 


Juif,{x,t),Ul{t) +Tv{t))dT v{t). 


We write Jx{f,ix, t) + sy{t) , U2{t)) as Jx and Ju{f{x,t),ui{t) + Tv{t)) as Ju- Then the differentiating 
||?/(t)|P with respect to t\ 

d 


(14) 


{Jx + Jx)dsj y{t) 

+ f [ Judr] y{t) + y^{t) ( [ Judr ) v{t) 


< y'^it) { Jx + Jxds y{t) + y'^{t)y{t) 


Judr ] v{t) 


Judr v{t) . 


+ 


Recall that XmaxiS) is the upperbound of the eigenvalues of the symmetric part of Jx over S, 
so Jx + Jx ^ 2Amax(*S')/. Therefore, (14) becomes: 


fL 

dt 


< {‘^Xmax{S) + l)||?/(t)|P -|- 


Judr v{t)\(^ 
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Let 2a = 2\max{S) + 1, M = max (|| t),tt)||), then equation (14) becomes 

u€U{I) 


|ll!/(*)P<2a||!,{()P + Af||o(''“2 
Integrating each side from tj_i to t where t < ti, we have: 


( 15 ) 


It follows that, 


||s(t)||" < j M\\v{T)fdT] ^ 


II'SWII < f ||»( 


T dT. 


(16) 


Using Lemma 4.5 to get the coefficients a[i] and M[i] in each time interval = 1... ,k, 

we will have: 

Theorem 4.6. The items in array a and M are a coefficients of the {K,hl(X))-local IS discrepancy 
function for the system (11). 

This theorem enables us to compute the {K,U{I))-local IS discrepancy function for each sub¬ 
system Ai. Although in the original definition we assume the IS discrepancy function is valid for 
any input signals ui,U 2 G U, in practice Ai can only take Aj’s outputs or states as inputs, which 
is bounded. Thus, [17] can still use {K,U{I))-local IS discrepancy function computed by this ap¬ 
proach. Furthermore, the {K,U(X))-local IS discrepancy function here can over-approximate the 
reachset of the systems in ( 11 ) with the input u being chosen nondeterministically in some compact 
set. 

5 Experimental Evaluation 

We have implemented the verification algorithm of Figure 1 and the ComputeLDF subroutine both 
with and without coordinate transformation in Matlab. The implementation and the examples are 
available from [12]. For simulation we use Matlab’s built-in ODE solver. The Jacobian matrix, 
an upper bound of the Lipschitz constant are given as inputs. In addition, the function to do the 
term-wise maximization of the error matrix is also given as inputs (see Section 3.2). We use the 
absolute error for ODE solver as the error bounds for simulation. The results presented here are 
based on experiments performed on an Intel Xeon V2 desktop computer. 

5.1 Comparison with other tools 

We compare the performance of our algorithm with two other tools, namely. Flow* [ 6 ] and HyCre- 
ate [22], for safety verification problem of nonlinear dynamical systems. We use seven benchmarks 
which are shown in Table 1 with time bound T = 10s. Flow* uses Taylor models for approximating 
reachtubes from a set of initial states. Currently, it returns “Safe” or “Unknown”, but not “Un¬ 
safe” . HyCreate uses the face-lifting approach of [ 8 ] and provides a intuitive interface for creating 
models. 
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Vanderpol, CoupledVanderpol, JetEngine, and Brusselator are commonly used, low-dimensional, 
nonlinear benchmarks. Sinusoidal tracking [21] is a 6 dimensional nonlinear designed as a frequency 
estimator. The Lorenz Attractor (row 7) is a well known chaotic dynamical system. Robot arm is 
a 4 dimensional nonlinear system described in [3]. The Helicopter is a high dimension linear model 
of a helicopter system from [13]. 

We have implemented verihcation algorithmwith and without coordinate transformation. Columns 
(^SimO) and (LDFO(s)) show the number of simulations and running time of our algorithm (Fig¬ 
ure 2) without coordinate transformation. In comparison. Columns (#Sim) and (LDF) are the 
results with coordinate transformation. Coordinate transformation provides tighter bounds, so the 
number of simulations and running time decrease under the same environment (i.e. same initial 
sets and unsafe sets). In row 10 and 11, we increase the time bound of the fixed-wing model to 
T = 50 and T = 100 respectively and the results show that the algorithm scales reasonably for 
longer time horizons. Flow* and HyCreate generate a single over-approximation of the reachtube 
from the initial set independent of the safety property. While our algorithm will rehne the initial 
sets when the reachtube intersects with the unsafe set. In all of these benchmarks, we make the 
unsafe set close to the reachtubes, to make the models safe yet it needs a lot of refinements to arrive 
at that conclusion. Overall, the proposed approach with coordinate transformation outperformed 
others in terms of the running time, especially in high dimensional benchmarks. The “N/A” in the 
table means the algorithm timed out at 30 minutes. Of course, our implementation requires the 
users to give the symbolic expression of the Jacobian matrix and term-wise maximization functions, 
while Flow* and HyCreate just needs the differential equations. Moreover, our implementation cur¬ 
rently handles only nonlinear dynamical systems, and both Flow* and HyCreate can handle hybrid 
systems. 

5.2 Properties of LDF 

We explore the behavior of the algorithm with respect to changes in the relative positions of the 
initial set and the unsafe set. We use the nonlinear model of the Robot arm system. We hx the 
point [1.5,1.5, 0,0] as the center of the initial set and T = 10 seconds as the time bound, and 
vary the diameter of the initial set (J) and the unsafe set (U : 0 > c), where 9 is the angle of the 
arm. The number of simulations used by the algorithm with coordinate transformation (^^^Sim), 
the diameter of the reach tube at the hnal time T (dia), and the total running time (RT) are shown 
in Table 2. 

From the first 5 rows in the Table, we see the expected behavior that for a fixed unsafe set, the 
diameter of the Reachtube decreases with decreasing 6 . This corresponds to the property that the 
discrepancy function /3{x,x',t) goes to 0 as the initial points x —>■ x', and therefore the error in the 
reachability computation decreases monotonically with the diameter of the initial set. Rows 4 and 
6-9 show that if we fix the size of the initial set, then as the unsafe set comes closer to the actual 
reachtube, the number of simulations increases and therefore the running time increases until the 
system becomes unsafe. As more refinements are made by the algorithm, the accuracy (measured 
by the diameter of the reachtube) improves. Similar trend is seen in rows 10-12, the algorithm will 
need more rehnements to find a counter example that shows unsafe behavior, if the unsafe set is 
close to the boundary of the reachtube. 

Next, we explore the behavior of the algorithm (with coordinate transformation) with large 
initial sets. We use the 7 dimensional model of a hxed-wing UAV. The initial sets are dehned as 
balls with different radii around a center point [30, 980,0,125,0,0, 30.4] and 6 in the first column is 
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example 

dim 

<5 

U 

#Sim 

LDF(s) 

#SimO 

LDFO(s) 

flow*(s) 

HyCreate(s) 

1 

Vanderpol 

2 

0.5 

x>2.0 

9 

0.378 

61 

2.01 

11.2 

2.776 

2 

Brusselator 

2 

0.5 

x>1.3 

21 

1.01 

85 

2.79 

11.8 

1.84 

3 

Jet Engine 

2 

0.4 

x>2.0 

5 

0.353 

61 

1.97 

8.74 

5.54 

4 

Robot arm 

4 

0.5 

x>2.5 

81 

4.66 

1159 

47.9 

169 

>300 

5 

CoupledVanderpol 

4 

0.5 

x>2.5 

41 

2.21 

1353 

54.2 

93 

49.8 

6 

Sinusoidal Tracking 

6 

0.5 

x>10 

185 

13.2 

753 

97.0 

258 

>300 

7 

Lorenz Attractor 

3 

0.02 

x>le4 

570 

13.99 

3105 

72.0 

53.4 

N/A 

8 

Fixed-wing UAV (T=10) 

7 

3 

x> 39 

321 

20.8 

N/A 

N/A 

N/A 

N/A 

9 

Helicopter 

28 

0.02 

x>4 

585 

67.7 

N/A 

N/A 

N/A 

N/A 

10 

Fixed-wing UAV (T=50) 

7 

3 

x> 39 

321 

99.8 

N/A 

N/A 

N/A 

N/A 

11 

Fixed-wing UAV (T=100) 

7 

3 

x> 39 

321 

196 

N/A 

N/A 

N/A 

N/A 


Table 1: Safety verification for benchmark examples, dim; dimension of the model; 6 : diameter of the initial set; U: unsafe 
set; #Sim: number of simulations with coordinate transformation; LDF: running time of our implementation (with coordinate 
transformation) in seconds; #SimO: number of simulations using algorithm in Figure 2; LDFO: running time of algorithm in 
Figure 2(without coordinate transformation)in seconds. 




















5 

U 

saftey 

#Sim 

dia 

RT(s) 

1 

0.6 

e >3 

safe 

17 

5.6e-3 

0.948 

2 

0.4 

e >3 

safe 

9 

3.6e-3 

0.598 

3 

0.3 

e >3 

safe 

9 

2.6e-3 

0.610 

4 

0.2 

e >3 

safe 

5 

1.8e-3 

0.444 

5 

0.1 

e >3 

safe 

1 

1.5e-3 

0.271 

6 

0.2 

e >2.5 

safe 

9 

1.7e-3 

0.609 

7 

0.2 

e >2.18 

safe 

17 

1.4e-3 

0.933 

8 

0.2 

e >2.17 

safe 

29 

l.Oe-3 

1.429 

9 

0.2 

e >2.15 

safe 

161 

9.2e-4 

6.705 

10 

0.2 

e >2.14 

unsafe 

45 

N/A 

1.997 

11 

0.2 

e >2.13 

unsafe 

35 

N/A 

1.625 

12 

0.2 

e >2.1 

unsafe 

1 

N/A 

0.267 


Table 2: Safety verification for a robot arm with different initial states and unsafe sets, safety: 
safety result returned by verification algorithm; 


the diameter of the initial sets. The unsafe set is defined as H > c, where H is the thrust of UAV. 
The time horizon is fixed at T = 10 seconds. As shown in Table 3, our algorithm can handle large 
initial set and high dimension systems. Although it may need many simulations (24001 covers), 
the algorithm terminates in 30 mins. All the results of this table are safe. 



(5 

U 

#Sim 

RT(s) 

1 

50 

H > 400 

24001 

1518 

2 

46 

H > 400 

6465 

415 

3 

40 

H > 400 

257 

16.33 

4 

36 

H > 400 

129 

8.27 

5 

20 

H > 400 

1 

0.237 


Table 3: Safety verification for a fixed-wing UAV with large initial sets. 


6 Related Work 

Simulation based verification has been studied in several papers recently [9, 1, 7, 18]. In [9] the 
authors introduce a general simulation based method for proving safety of arbitrary continuous 
systems. The novelty of their approach consist in the use of sensitivity analysis, where the sensitivity 
matrix with respect to initial state xq at time t is defined as s^o — ■ It is shown that 

Sxo(t) = Jf{xo,t)sxg(t) and Sx^it) can be solved by efficient solvers. Then ||sxo(^)ll<^ is used to 
bound the distance ||^(x,f) — C(xo,t)|| for x E ^^(xo) at time t. It is shown that this upperbound 
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holds for linear time varying systems. For general nonlinear systems, ||s2;Q(t)||(5 has a quadratic 
error term with respect to 6 that requires further analysis. Thus, this technique is sound for linear 
system but does not provide any formal guarantees for nonlinear systems ([9], page 13). In [7] this 
technique is extended to nonlinear systems subject to disturbances as inputs and uncertainty in the 
initial conditions to obtained an approximation that ignores the higher order terms. In contrast, 
in Section 3 and Section 4 we have provided a strict over-approximation of Lipschitz continuous 
systems with respect to uncertainty in the initial conditions and uncertainty in the input signals. 
In [18], the authors provide several approaches to capture the upperbound of the distance between 
two trajectories for linear systems and some polynomial systems. 

In [15] the authors present a convenient implementation of sensitivity analysis in the Simulink 
software. Again, the trajectory sensitivity matrix can only be used as a linear approximation for 
a perturbed trajectory , instead of over-approximation of the reachset. In [1] the authors provide 
a different approach by linearizing the nonlinear system locally, and bounding the linearization 
error by Lagrange remainders. The original definition of discrepancy function can be seen as a 
generalization of the incremental stability [2]. The incremental Lyapunov function can be used as 
discrepancy function when a system is incrementally stable. An incremental Lyapunov function- 
based approach is used in [14]. Here the authors go much further and construct a finite symbolic 
model that is approximately bisimilar to the original switched system. Our approach bypasses the 
incremental stability requirement by focusing on bounded time analysis. 

Contraction in [19] is defined as the region in which the eigenvalues of the symmetric part of 
the Jacobian is uniformly negative. The authors use “virtual displacement” to get the result, while 
we get the upperbound of the eigenvalues of the symmetric part of the Jacobian directly from the 
generalized mean value function. Contraction metrics introduced in [19] is also used in [10] to 
perform sound and relative complete analysis of nonlinear systems. 

7 Conclusions and Future Work 

In this paper, we present an algorithm ComputeLDF to compute local discrepancy functions, which 
is an upperbound of the distance between trajectories starting from an initial set. The algorithm 
computes the rate of trajectory convergence or divergence for small time intervals and gives the rate 
as coefficients of a continuous piece-wise exponential function. The local discrepancy we compute 
satisfies the definition of discrepancy function, so the verification algorithm using ComputeLDF 
as a subroutine is sound and relatively complete. We also provide a coordinate transformation 
method to improve the estimation of rates. Furthermore, we extend the algorithm to compute 
input-to-state discrepancy functions. 

In the future, we plan on using more rigorous ODE solvers like [5] and embedding the algorithm 
in verification tools like C2E2 [10] for safety verification of hybrid systems. 
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A Appendix: Proofs of Lemmas 


Proof of Lemma 3.3: 

In this proof, the i’s in subscript correspond the the components of the vector functions. 
For any t E [0,1], i E {1,... , n}, we define gi{t) := fi{x + tr). Then we have 


fi{x + r)-fi{x) = gi{l) - gi{0) = [ 

Jo 


^ dgiit) 


dt 


dt. 


(17) 


Using the chain ruie of gradient, we have 

dgiit) 


dt 




dix + tr) 


u=x-\-tr 


dt 


= ^fiiu)\u=x+tr-r = Y1 


dfiiu) 


i=i 


dui 


'31 


(18) 


u=x-\-tr 


where Vfi{u) = ,..., is the gradient of function /j. Substituting (18) in (17), 

we have: 


fiix + r) - fi{x) 



dfiju) 


u=x-\-sr 


ds 



ds 


u=x-\-sr 


Since Jfix + sr) is the matrix consisting of the components of -4^ j the iemma hoids. 

d u=x-\-sr 

Proof of Theorem 3.4. 

This theorem is estabiished by the minimax characterization of the eigenvaiues. Let A = A + E, 
and iet Aj(A), XiiE), Ai(A) denote the eigenvaiues of A, E and A respectiveiy , where aii three sets 
are arranged in non-increasing order. By the maxmin therorem we have 

Xk(A) = min ( max pziv) ) 

^ ^ dimV=n-fc-|-l 

Which can aiso be written as 


Xk{A) = min max(x'^Ax) 
x^x = 1,pjx = Oii = 1,2,3,... ,k — 1) 

Hence, if we take any particuiar set of pi, we have for aii corresponding x, 

Afe(^) < max(x"^Ax) = max(x'^Ax -|- x'^Ex). 

If U'^AU = A = diag(Aj(A)) and U is the orthogonal matrix, then if we take pi 
relations to be satisfied are 

0 = pf x = ef yii = 1,2 ,... , k - 1) 


(19) 
Uci the 
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With this choice of the pi then the first k — 1 components of y are zero, and from equation (19) 
we have 

n 

< max(x'^^x + Ex) < max(^^ \i{A)y‘f + x^ Ex) (20) 

i=k 

However, 

n 

j;\(H)y2<A,(H) (21) 

i=k 

while 

x^Ex < Ai(i?) (22) 

for any x. Hence the expression in brackets of equation (20) is not greater than Afc(H)+Ai(i?) for any 
x corresponding to this choice of the pi. Therefore its maximum is not greater than Afc(H) + Xi{E) 
and we have 

Afc(i)<Afc(H) + Ai(S) (23) 

Since A = A+{—E) and the eigenvalues of —E in non-increasing order are —Xn{E), —Xn-i{E), ..., — Ai(£'), 
and application of the result we have just proved gives 


Afc(j4) < Xk{A) (—Xn{E)) or Xk{A) > Xk{A) + Xn{E) (24) 


Thus we have 

Afc(^) + Xn{E) < Xk{A + E) < Xk{A) + Xi{E) 

The relations (23) and (24) imply that when E is added to A all of its eigenvalues are changed 
by an amount which lies between the smallest and greatest of the eigenvalues of E. Note that we 
are not concerned here specifically with small perturbations and the results are not affected by 
multiplicities in the eigenvalues of A, E and A + E. 
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